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Introduction 

This report summarizes progress made from 1 June 1986 through 30 
November 1986 in developing an effective algorithm for estimating vegeta- 
tion cover at subpixel resolution. An important application is the 
estimation of soil hydraulic properties in natural water-limited vegeta- 
tion systems using Eagleson' s "climatic-climax" hypothesis (Eagleson, 
1982, Eagleson and Tellus, 1982). It is also critical to the 
subgrid-scale parameterization of heat and moisture fluxes across the 
atmosphere-landsurface interface in the development of global-scale 
numerical atmospheric general circulation models (GCMs). Such models are 
essential for the study of climate change and for large-scale 
environmental impact assessment. Calculation of the fluxes requires 
knowledge of effective areal averages of soil hydraulic properties and of 
vegetation characteristics at the mesoscale. 

An earlier report (Jasinski and Eagleson, April 1986) concluded that 
only a few methods for resolving mixed pixels were available, and they 
had not been tested for the special case of determining percent canopy 
cover. Two formulations were proposed: One using the normalized vegeta- 

tion index, a second involving a simplified expression of the radiative 
transfer equation. 

That report also includes details of the database acquired for the 
project study area near Taos, New Mexico. Landsat MSS data were obtained 
from the Bureau of Land Management (BLM) Branch of Remote Sensing, 

Denver, Colorado. Ground truth in the form of 1:3000 color aerial photo- 
graphs were borrowed from the BLM Taos Resource Area Office, Taos, New 

» 

Mexico. ' 
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The present report summarizes the various approaches relevant to 
estimating canopy cover at subpixel resolution. The approaches are based 
on physical models of radiative transfer in non-homogeneous canopies and 
on empirical methods. The effects of vegetation shadows and topography 
are examined. Simple versions of the model are tested, using the Taos 
Study Area database. Emphasis has been placed on using relatively simple 
models requiring only one or two bands. Although most methods require 
some degree of ground truth, a two-band method is investigated whereby 
the percent cover can be estimated without ground truth by examining the 
limits of the data space. Future work is proposed which will incorporate 
additional surface parameters into the canopy cover algorithm, such as 
topography, leaf area, or shadows. The method involves deriving a 
probability density function for the percent canopy cover based on the 
joint probability density function of the observed radiances. 

Homogeneous Canopy Reflectance 

The radiation reflected from horizontally homogeneous canopies 
results from the scattering and reflectance properties of the plant 
components and soil background. These properties are both geometric and 
biophysical in nature and thus depend on the species, maturity and over- 
all health of the plant. Geometric plant properties include plant archi- 
tecture, total leaf area, and leaf orientation and distribution. Bio- 
physical properties allocate the radiative energy absorbed by the plant 

i 

to important metabolic processes including photosynthesis, respiration, 
and transpiration. Those biophysical properties are manifested in terms 
of leaf color, transparency, temperature, and shape and orientation. 

Plant properties can vary daily and seasonally, in response to soil 
moisture and nutrients, and to meteorologic and climatic conditions. 
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Numerous radiative transfer models for horizontally homogeneous 
canopies have been developed in terms of various plant properties and 
background soil reflectance. Typically, homogeneous canopies have been 
modeled as a diffusing medium with absorbing and scattering properties. 
Excellent reviews of these models are provided by Smith (1983) and Ross 
(1984). Suits (1972) envisioned a plant canopy as an infinitely extended 
plane-parallel medium with homogeneous geometric properties. Verhoeff 
and Bunnik (1981) extended the soils model to include the effect of leaf 
angle distribution. Dickinson (1983) applied the two-stream approxima- 
tion for radiation transfer in the atmosphere (Meador et al., 1980) to 
plant canopies employing the leaf area index (LAI) as a measure of the 
optical depth. Recent literature has increased the sophistication of 
those earlier models to include the modeling of bidirectional reflectance 
(Walthall et al., 1985, Chen, 1985, Simmer et al., 1985, Vanderbilt et 
al., 1985, Strebel et al., 1985, Gerstl et al., 1985, and Reyna et al., 
1985). Attention has also been focused on the invertibility of canopy 
reflectance models for estimating plant parameters such as LAI, biomass, 
and leaf angle distribution (LAD) (Goel et al., 1983, Goel et al., 1984a, 
Goel et al., 1984b, 1984c, 1984d, Lang et al., 1985). Goel et al. (1984) 
have shown that such parameters can, in principle, be estimated using 
only canopy reflectance measurements at several viewing angles. 

Because such theoretical models are often cumbersome to use, semi- 
empirical formulas for the total radiation fluxes have been proposed for 
practical applications. The attenuation of radiation as it passes 
through a plant stand has been typically described in terms of some form 

of Bouguer's Law such as that proposed by Monsi and Saeki (1953), or 

s' 

I(L) = I exp(-yL) 
o 


( 1 ) 
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where 

I = intensity of radiation at top of canopy 

I(L) = intensity of radiation at a penetration level 
associated with leaf area index, L 
y = experimentally determined extinction coefficient 
The plant reflectance, r y , is thus 


r 


v 



( 2 ) 


Other formulas which account for the different attenuation of 
penetrated and scattered radiation, as well as absorbed radiation, have 
been proposed (Ross, 1981). 

Actual reflectances of plants in the visible and near infrared are 
well documented, especially for crops (Smith, 1983, Myers, 1983). Canopy 
reflectance is highly wavelength dependent as shown in Exhibit 1. It is 
typically low in the visible spectrum (< 30%) and higher in the near 
infrared region (> 50%). Canopy reflectances are generally much lower 
than those measured for individual leaves (Dickinson, 1983). Kondratyev 
(1969), Gates (1980), Ross (1981) and Iqbal (1983) provide summaries of 
reflectances for various natural vegetation covers. 

Non-Homogeneous Canopy Reflectance 
Non-homogeneous canopies are those which contain several vegetation 
or soil types within the level of resolution being investigated. Most 
natural landscapes will vary both horizontally and vertically in species 
and/or vegetation density. Modeling this situation has received 
considerably less attention than have homogeneous canopies. Statistical 
techniques have been employed for classifying landscapes. However, such 
methods require the identification of training samples and therefore can ^ 
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not be adopted for natural landscapes in which all target pixels possess 
unique spectral characteristics. 

Radiative transfer models for non-homogeneous canopies have been 
developed by extending homogeneous canopy models and including 
three-dimensional scattering functions (Ross, 1981, Kimes, 1984). Such 
models can be solved in a few cases where plant distribution is periodic, 
such as for sown crops. Applications to natural systems is impractical 
if not impossible as it would require knowing a priori an inordinate 
number of plant parameters. 

The problem can be somewhat simplified by recognizing that for most 
hydrologic or climate modeling purposes, it is generally sufficient to 
know only the outgoing radiance at the surface of a canopy, and not the 
complex scattering and absorption phenomena within the canopy. 

Researchers have investigated mixed pixels by examining only the 
surface reflectivities of the individual components. Two such approaches 
are summarized below. 

Proportion Estimation 

A simple approach has been to assume that the surface reflectivities 

of individual plant clusters and soils within a pixel are constant, and 

that the total spectral response is a linear combination of the 

individual spectral responses of its components (Horwitz et al., 1971, 

Nalepka et al., 1972, Work, 1974, McCloy, 1980, Dozier, 1981, Ungar et 

al., 1981, Chhikara, 1984). The total radiance emitted from a pixel, I , 

X 

containing n cover types can be expressed 


n 


I 


i«l 


m i I Xi 


(3) 
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where 

= fraction of area covered by cover type i; 

I * radiance emitted from cover type i in band X • 

A -L 

For the simple case of vegetation and soil cover, each with constant 
reflectivities, Equation (3) becomes: 

\ • mI Xv + <4) 

m = percent vegetation cover 

I, ,1 * radiances emitted from vegetation and soil, 

Xv’ xs 6 

respectively, in band 

Perhaps the earliest development of proportion estimation can be 
attributed to Horwitz et al. (1971) and Nalepko et al. (1977) who also 
termed the method "mixtures estimation". One of its first applications 
was in identifying subpixel scale ponding and wet marshes in glaciated 
prairies (Work, 1974). 

McCloy (1980) later proposed that under conditions of negligible 
canopy transmission or multiple reflection, the response proportions of 
the various land covers will closely approximate the physical proportions 
of each type of cover. He suggested that up to four sub-pixel categories 
be used including three levels of vegetation greenness cover and one soil 
background cover. Ungar et al. (1981) reported limited success 
delineating forest canopy types in Maine using a similar approach which 
they termed the "Fanning algorithm". The functional area was determined 
by minimizing the error between observed and theoretical radiances. 

Dozier (1981) also proposed such a method using two infrared bands 
for differentiating radiant temperature fields of sub-pixel spatial 
resolution. Corrections for atmospheric effects were included. 
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Geometric Models 

In an extension to the linear proportion estimation model described 
above, some investigators have considered the shadow cast by vegetation 
as an additional component to the total radiance. These models abstract 
clumps of vegetation as three-dimensional geometric shapes on horizontal 
surfaces with constant reflectivities (Otterman 1981, 1984, Otterman et 
al., 1984, Strahler et al., 1981, Li et al., 1985). 

Otterman (1981, 1984) and Otterman et al. (1984) envisioned forests 
and desert vegetation as thin vertical cylinders. While accounting for 
the shadowing effects of vegetation clumps, the formulation assumes that 
reflection from the top of the clumps is negligible. This model is thus 
not directly applicable to the determination of percent vegetation cover. 

Strahler et al. (1981) and Li et al. (1985) modeled conifer forests 
as randomly located cones of similar shape and random height. They 
determined from simple geometry the shadow cast by the cones on the soil 
background or on other cones. The total radiance emitted by a pixel was 
assumed to consist of four components: illuminated background, 

illuminated cones, shadowed background and shadowed cones. Vegetation 
parameters including percent cover and average tree height were then 
estimated using assumed values of component reflectivities. 

A simple version of the geometric model is presented below. 

Assuming the radiance emitted from shadowed areas negligible, the total 
radiance from sparse conifer trees, I , modeled as cones on a horizontal 

A 

soil background with constant reflectivities is. 


.W r i» ( ? + ^ ' Vt + 
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0 - solar zenith angle 

a = ratio of tree height to base radius, 
assumed to be constant 

I = incident radiation in band 

Xo 

The model can be extended to include non-horizontal surfaces and 
overlapping of shadows# Li et al. (1985) further assumed a poisson 
distribution of the cones and random overlap of shadows. 

Satellite Observed Radiance 

In order to derive canopy parameters from satellite data, it is 
necessary to couple the radiative transfer through the atmosphere to that 
through the plant canopy. As a starting point, it is useful to review 
the atmospheric models which have been developed for the visible and near 
infrared regions. 

The theoretical nadir radiance observed by a satellite has been 
derived by several authors (Dave 1980, Otterman et al., 1979, 1980, and 
Otterman, 1978, 1981) for a pure pixel in a background of different 
reflectivity. The nadir radiance is expressed in terms of three 
components. The direct beam from the object pixel which has been 
attenuated due to atmospheric effects, , is 

L nrx * rG x <6) 

where r ■ reflectivity of object pixel 

G “ global surface irradiance on object pixel 

A 

(Otterman, 1978) 

t = total vertical optical thickness, a function of \ * 

A 

The total optical thickness is a sum of optical' thicknesses due to 
Rayleigh scattering by gas molecules, x R * and Mie scattering by aerosols, 
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Tw f also functions of 
M 

Assuming a Lambertian surface, the portion of diffuse radiance from 
the surrounding vicinity which is scattered to the satellite by the 
atmospheric column above the object pixel is 

tt/2 

L naX = (aG^/TTi x )f [l-expC-x^/cos <|>)]cos 4 > [ t aR p xR (<I > ) +t xM P Xm(^ > ^ 2ir sin 4 ,d <l > 
0 (7) 

where 

a • average reflectivity of the surrounding vicinity 

P „,P „ = phase functions associated with Rayleigh and 
XR AM 

aerosol (Mie) scattering, respectively 
= zenith reflection angle. 

Finally, the radiance scattered from the direct solar beam back to the 
satellite is written 

p [ 1-exp {-T^(l+sec 0 )} ] [ t xr P xr(^0° “9 ) + (8) 

L ndX " t x (1 + v ) 

where p = cos 0 

0 = solar zenith angle 

The total nadir radiance is the sum of equations (6), (7) and (8). The 
amount of area to be considered in determining the average surrounding 
reflectivity is discussed by Otterman et al. (1979). Otterman et al. 
(1980) reduced the above equations to a simpler form in the case of an 
optically thin atmosphere. 

Subpixel Resolution of Satellite Observed Radiance 
The atmospheric model is coupled to the canopy models through the 
reflectance parameter, r, in the direct beam equation. For homogeneous 
canopies, the direct beam equation can be rewritten for example, 
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L 


nr A 




(9) 


For non-homogeneous landscapes, the proportion estimation equation (3) 
can be included into the direct beam equation yielding 


where 


L nrx’T eXp W K Vxil 


r^ ■ reflectivity of component i 

* fraction of pixels exhibiting reflectance r^. 


( 10 ) 


Finally, the conical geometric model incorporated into the direct beam 
equation yields 


nr A 


For clear skies in which the atmosphere is considered optically 

thin, direct beam radiation can constitute 80 - 90% of the total 

radiation observed by a satellite. The more complex scattering terms can 

often be neglected, which leaves a relatively simple expression of 

satellite observed radiance in terms of G, , x, , m, r, , r and v. A 

X A Av’ Av 1 

'straightforward application of proportion estimation for a two component 
system, i.e., soil and vegetation, would include equation (10) with four 
unknowns; m, r r , and the term G exp(-x.)/ir. 

A v A S A A 

Inclusion of shadow effects from a conifer forest (modeled as cones) 
would involve an additional unknown, y. 

The application of equations (10) or (11) to the determination of m 
depends on the data set being used and on the number of a priori assump- 
tions one is willing to accept about other parameters, such as ^ 

rXv, rXs, and T X. The use of multispectral bands can help reduce the 
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number of unknowns, especially if reasonable assumptions are made 
regarding the relationship between a given parameter in different bands. 
The results of several variations of these models tested with field data 
are provided in a later section of this report. 


Two Band Method 

The ideal canopy cover algorithm would not only provide percent 
coverage of individual pixels, as well as several parameters such as 
vegetation and soil reflectivities, but would do so without the need for 
ground truth. The following procedure was investigated as a step toward 
achieving that goal. 

It has been observed in a manner similar to previous authors (e.g., 
Kauth and Thomas, 1978), that MSS data plot in characteristic shapes in 
the two, three, or four-band space. For the MSS bands 4 and 2, with band 
4 plotted along the ordinate and band 2 along the abscissa, the data 
typically fall in a triangular-like shape in the upper half plane with 
curved sides and a flat base and with sometimes a "tasseled cap". This 
has indeed been true for the Taos Study Area whether one plots the entire 
region (all types of vegetation) as shown in Exhibit 2, or segments of 
the region, as shown in Exhibit 3 for the valley area, and Exhibit 4 for 
the foothills region. 

It can be shown using the direct beam equation that such character- 
istic triangular shapes can be explained in terms of the percent coverage 
of vegetation as well as the absolute values of the vegetation and soil 
reflectances. The starting point is the expression for direct beam, 
equation 10, written in terms of reflectances for MSS bands 2 and 4, or 


“ (r 2v- 


r 2s> + r 2s 


( 12 ) 
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r. = m(r. - r. ) + r. 

4 4v 4s 4s 

where 

r 2v* r 4v = vegetation reflectance in MSS bands 
2 and 4, respectively, 

r„ ,r. = soil reflectance in MSS bands 

2s 4s 

2 and 4, respectively, and 

r^ ,r^ = total observed reflectance of pixel 

in bands 2 and 4, respectively. 


(13) 


A simple triangular shape can be achieved by setting r 2 v > an< * r ^ v 
constant, and letting r^ and r^ g vary over a range of values. A line 
representing m = 0, or "soil line" (or "background line”) will form the 
base of the triangle stemming from the origin. Lines for values of m > 0 
will plot above but parallel to the soil line. The value for m=l, or 
full coverage, will plot as a single point at the apex of the triangle. 
For instance, a simple triangle is illustrated in Exhibit 5 for the case 
of a landscape with a single vegetation species with constant 
reflectivity. That shape was drawn using « 0.10, r^ v = 0.40, and 

soil reflectivities falling in the intervals 0.08 < r£ g j< 0.40 and 0.10 j< 
r^ g 0.50. The only constraint is that r 2 s ^ v i is rema ^ n approximately 
constant which is a reasonable assumption. 

The effect of changing the reflectivity values is shown in Exhibit 6 
(changes of vegetation reflectivity) and Exhibit 7 (changes in soil 
reflectivity). For instance. Exhibit 6 indicates that increasing r^ v 
will tend to flatten the triangular data space by moving the top of the 
triangle toward the soil line. An increase in will cause the 

opposite effect. A change in either r^ v or r^_ will have no effect on ^ 
the soil line as one expects. 
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The simple model yielding a triangular shape can be extended to 
include the effects of varying vegetation reflectivity and shadows. For 
instance, by assuming that vegetation reflectivity is linearly related to 
percent cover, the triangular shape takes on curved sides as shown in 
Exhibit 8* The inclusion of shadows caused by conical figures will cause 
the triangle to take on a tasseled cap as shown in Exhibit 9. While 
other variations in plant parameters might cause similar effects, the 
above examples illustrate how the characteristic shape of the data space 
can be explained in terms of physical attributes of the components of the 
pixel. Application of the two-band model to field data is provided later 
in this report. 

Empirical Estimators of Vegetation Parameters 
Numerous vegetation indices have been proposed in recent years as 
qualitative indicators of green vegetation. The purpose has been to 
reduce the several multispectral bands to one value to estimate vegeta- 
tion parameters such as biomass, leaf area index, or percent cover. 

Perry et al. (1984) summarize the many different vegetation indices and 
describe their relationship to each other. Three such indices are the 
normalized vegetation index, the perpendicular vegetation index, and the 
Kauth-Thomas Greenness index. 

Normalized Vegetation Index 

Of the many indices proposed, the normalized vegetation index (NVI) 
has evolved as a practical popular choice for use in regression with 
vegetation parameters. It is of the form 
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NVI 


r iR ~ r VIS 
r lR + r VIS 


(14) 


where and r^g represent the pixel reflectivity in the IR and visible 
ranges, respectively. Low NVI indicates low vegetation amount, whereas 
high NVI indicates either high vegetation amount or high productivity 
(Curran, 1980). Sellers (1985) discussed the functional relationship 
between the normalized VI and several vegetation parameters, including 
the leaf area index. Tucker et al. (1983) correlated the NVI (computed 
from NOAA's AVHRR data) to actual biomass obtained from field sampling in 
a semi-arid region of Senegal, West Africa. The effect of soil back- 
ground on NVIs computed from hand-held or airborne radiometer data was 
investigated by Elvidge et al., 1985, and Huete et al., 1985. 

Physical Basis. While the behavior of NVI can be qualitatively 
described in terms of chlorophyll absorption in the red region and 
scattering due to leaf area index in the near infrared, it can also be 
explained in terms of percent cover using the direct beam equation 
written in terms of reflectivities (12) and (13). Inserting (12) and 
(13) into (14) yields, 


NVI = 


ml<r 4v~ r 4s> - (r 2v~ r 2s )] * r 4s ~ r 2s 
“ I(r 4»- r 4s > + <r 2»- r 2s )] + r 4s + r 2s 


(15) 


It is now noted that theoretically, the reflectance of any pixel 
composed of only one type of soil and vegetation must fall between the 
limits prescribed by the pure canopy and pure soil reflectivities shown 
in Exhibit 1. Using typical values for r^ g , r 2 s » r ^ v and r 2 v as indica- 
ted, it is reasoned that the expression (r^- r^ g ) generally a posi- 
tive value, whereas (r„ v - r 2 s ^ generally negative. As m increases. 
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the numerator becomes more positive, while the denominator decreases. 

Thus NV1 should increase with increasing percent vegetation cover, as has 
been observed. 

Perpendicular Vegetation Index 

Richardson and Wiegard (1977) proposed the perpendicular vegetation 
index (PVI) as a measure of plant development. Application of this index 
first requires the establishment of a background soil line by linear 
regression of MSS bands 2 and 4 using bare soil pixels. The soil line is 
thus a straight line stemming from near the origin. The PVI is the 
perpendicular distance from the soil line to the actual data point which 
contains vegetation, and is defined, 

PVI - l(r g2 - r p2 >2 + (r g4 - V 2, ‘ /2 (16) 

where 


r g2’ r g4 


r p2’ r P 4 


reflectances of soil background in bands 
2 and 4, respectively, corresponding to 
the data point. 

reflectances of data point in bands 2 and 4, respec- 
tively, perpendicular to r ^ and r ^ on the soil 


line. 

Richardson et al. (1977) regressed PVI with percent cover of sorghum 
with a correlation coefficient of 0.57. Theis et al. (1984) studied the 
effect of vegetation and soil moisture on PVI. Rosenthal et al. (1985) 
recently used the PVI to investigate crop height and biomass. 

Physical Basis. The behavior of the PVI can be readily demonstrated 
in terms of percent cover using the two-band method presented earlier, 
although the methods are slightly different. By setting r^ v and r^ v con- 
stant and letting r^ g and r^ g range from dark to bright values, a tri- 
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angular data space can be drawn. The base of the triangle, or soil line, 
stems from the origin and represents the condition m • 0. As m increa- 
ses, the data space moves perpendicularly from the soil line. 

Thus, the PVI performs remarkably well in describing the physical shape 
of the data space when explained in terms of percent vegetation cover. 
Kauth-Thomas Greenness Index 

Kauth and Thomas (1926) applied Gram-Schmidt orthogonalization to 
the original four Landsat bands resulting in a new four-dimensional space 
termed "tasseled cap”. The procedure, which is similar to principle 
components except in the order of calculations, essentially rotates the 
data so that most of the variability can be explained in terms of four 
indices; greenness (GI) , brightness (BI), yellowness (YI), and nonsuch 
(NI). The first two of these indices are defined, 

BI - 0.332 DN1 + 0.603 DN2 + 0.676DN3 + 0.263 DN4 (17) 

GI - -0.283 DN1 - 0.660 DN2 + 0.577DN3 + 0.388 DN4 (18) 

where DN1, DN2, DN3, and DN4 represent the digital counts of the three 
visible and one near infrared bands of the MSS scanner, respectively. A 
similar set of equations has been developed for the Thematic Mapper 
(Crist, 1983, Crist et al., 1984). 

The Kauth-Thomas Transformation has been used by numerous investiga- 
tors to model various crop parameters including crop development, mois- 
ture stress, yield and crop classification (Ezra et al., 1984). Huete et 
al. (1984, 1985) in a series of small scale experiments of wooden boxes 
filled with soil, showed high correlation of GI with percent cover. 

Musick (1984), however, using Landsat MSS data over New Mexico, was 
unable to achieve consistent differentiation between arid rangeland covert' 
changes using GI. 
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Application of Canopy Cover Algorithms to Taos Study Area 

Several different canopy algorithms described above have been tested 
at a preliminary level for a site located near Taos, New Mexico. The 
database consisted of Landsat MSS data, and 1:3000 aerial photographs, 
both supplied by the Bureau of Land Management. Details of the test site 
selection, and data acquisition and processing are provided in Jasinski 
and Eagleson (1986). 

The Taos Study Area is outlined in Exhibit 10. The land includes a 
wide variation in surface relief, ranging from flat plains to rolling 
foothills, to detached high ridges. Elevation ranges from 6,000 to 
10,000 feet. Vegetation tends to follow the topography. The lower flats 
are covered with blue grama and wheatgrass grasslands, and snakeweed, 
rabbitbrush and sagebrush shrublands. Pinyon- j uniper woodlands are found 
in the rolling foothills. At the higher elevations, there is ponderosa 
pine, spruce, fir and aspen. Percent cover ranges from nearly 0 to 100 
percent, with the majority of the area 40 to 60 percent covered. At 
least two trends in percent vegetation cover can be readily observed. 

They are, first, a decreased percent vegetation cover with decreasing 
altitude, and second, a less dense cover on south-facing slopes compared 
to north facing slopes at the same elevation. 

Regression with Normalized Vegetation Index 

Three variations of this approach were tested by regression with 
ground truth obtained from aerial photographs. A total of 116 pixels 
were used. The first variation involves using the normalized vegetation 
index with reflectivities as defined in equations (12) and (13). A 
second variation uses the normalized vegetation index defined in terms of 
actual integer DN values instead of reflectivities. The third approach 
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uses actual radiances. Details of these earlier studies were reported by 
Jasinski and Eagleson (1986), Only the results are provided in Table 1 
below and shown on Exhibits 11, 12, and 13. 


Table 1 


Normalized Vegetation Index 


versus Percent Cover 


NVI Variation 


m 


R 2 


VI = 


V r 2 


r V r 2 


m - 2.36 VI ' + 6.68 
r 


0.56 


VI 


2xDN. - DN„ 
4 2 

DN - 2xDN.+ DN. 

4 2 


m = 1.99 VI + 0.95 

DN 


0.61 


R 4" R 2 

VI R " * 100 


m = 3.06 VI - 111 


0.58 


The results indicate that about 60% of the change in NVI can be 
explained in terms of percent vegetation cover, no matter which variation 
is used. The method is quick, requiring little analysis or computational 
time. 

There are several limits to the use of NVI for estimating percent 
cover. First, the procedure needs to be calibrated with ground truth. 
Second, its accuracy is dependent on the vegetation and soil reflectivi- 
ties in each band being relatively constant, and also limited to r^ g > 
r^ v and r^ v > r ^ s * For regions which are relatively uniform in 
vegetation type, however, the method provides good results of percent 


cover. 
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Regression Using Direct Beam Equation 

Assuming that the landscape consists of only two cover types, soil 
and vegetation, equation (10) can be rewritten, 


m 


1; 


nr\ 


(r -r ) l G 

XV XS -Aexp(-. ) 
TT A 


‘"As] 


(19) 


For constant reflectivities, m is approximately linearly related to the 
observed radiance. With this in mind several linear regressions were 
carried out with m as the dependent variable and the observed radiances 
as dependent variables. The results are shown in Exhibits 14 and 15 and 
summarized in Table 2. 

Table 2 


Summary of Linear Regressions 


Regression Equation R 

m = 126 - 199R2 0.53 

m = 377 - 138R4 0.10 

m = 143 - 189R2/ cos 0 0.42 

m = 120 - 197 R2 cos & 0.53 

m = 135 - 185 R2 cos g/cos 0 0.44 


where R2 and R4 are the observed radiances in MSS bands 2 and 4, respec- 
tively. 

Two modifications to the above linear regressions were made. First, 
since is a function of the cosine of the zenith angle 0, equation 
(19) can be rewritten 


m 



V 


L , 
r nrA 

LG' COS0 

A 


ex p (-Tx) 



(20) 


IT 
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t 


Second, slight corrections to the observed radiance can be made in the 
case of non-horizontal surfaces. For surfaces with average slope g, the 
equation for m can be rewritten 


m 


G T 


L % cosg 
nr\ p 


(r -r ) L G T exp(-T ) xsj 

Xv xs X cose ^ X 


( 21 ) 


For the above two cases, ground slope and azimuth of individual pixels 
were measured from USGS 7.5-minute topographic maps. Zenith angle was 
computed using the following formula (Iqbal, 1983), 


Cos 0 * (sin^ cosg - cos<|> sing cosy) sing 

+ (cos<j> cosg + sin<j) sing cosy) cosfi cosoj 

+ cos5 sing siny sinu> (22) 

6 = solar declination 

a) = hour angle 

y = surface azimuth angle 

g = average slope of pixel 

Solar declination and hour angle were estimated from knowing time of 
Landsat overpass. 

The results of the linear regressions including the new parameters 
are surprisingly poor, as indicated on Exhibits 16, 17 and 18, and on 
Table 2. In no case is the correlation improved. Part of the explana- 
tion for the poor correlation may simply be due to the bidirectional 
reflection characteristics of the soil. A likely explanation may also 
simply be the inaccuracies introduced by measuring small distances off 
the topographic maps. At the 1:24000 scale, pixels are less than 0.2 cm? 
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in area and small inaccuracies in measurement or pixel registration can 

cause serious error in the regression analysis. 

Use of the direct beam equation for estimating m has required a set 

of ground truth. Having the ground truth, it is theoretically possible 

to estimate the quantities ir/[(r, - r, )G, _ exp(—r )] and r /(r -r ) 

M Xv As At * A As Av Xs' 

from the coefficients of the linear regression. If functional relation- 
ships can be established between the individual parameters for each band 
(e.g., r 2 ^“ f(r^)) it is a ^ so possible, at least in theory, to calculate 
all the reflectivity values r 2 v > r 4 v > r 2 s > r 4 s as we ii as the quantities 
G 2 t exp(~T 2 ) and G^ t exp(-x^) using data from bands 2 and 4. 

Regression with Kauth-Thomas Indices 

The Kauth-Thomas greenness and brightness indices were computed 
and then regressed with actual percent cover obtained from the aerial 
photographs. The results which are somewhat poorer than anticipated, are 
shown in Table 3. They indicate, contrary to expectation, that bright- 
ness appears to explain more of the variation in m than greenness. 

Table 3 

Kauth-Thomas Indices 


Index 


R2_ 

Greenness , 

GU 

0.24 

Brightness, 

BI 

0.39 

The regression analys 

es yielded. 


m = 10.47 

+ 0.14 GI 



m = 82.82 - 0.48 BI 
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Multiple Linear Regression 

Multiple linear regressions were carried out with the same data set 
as in previous cases with percent cover as the dependent variable and the 
MSS band observations as independent variables. The two cases examined 
were m vs, DN1 and DN2, and m vs. DN1, DN2, DN3 and DN4. Once regression 
coefficients were obtained, theoretical percent cover obtained from the 
multiple linear regression analysis was regressed with actual percent 
cover in order to compare correlation coefficients with other methods. 

The results for the second case are shown on Exhibit 19 (for regression 
with four bands) and summarized on Table 4. 

Table 4 

Results of Multiple Linear Regression 

Regression Equation R 2 

m = -2.25 DN2 + 0.70 DN4 + 74.97 0.53 

m - -2.07 DN1 - 0.62 DN2 + 0.20 DN3 + 0.63 DN4 + 72.25 0.58 

As expected, there is negative correlation with the visible bands 
and positive with the near infrared. It is also noted that the addition 
of bands 1 and 3 only contributes an increase of 0.05 in R 2 . 

Two-band Model 

The two-band model as presented was applied to the foothills segment 
of the scene where the type of vegetation was believed to be relatively 
homogeneous. Two cases were investigated, those assuming constant and 

s' 

variable (with m) vegetation reflectivities. 
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Application of the model consists of the following steps: 

1) All data points within the segment are plotted. 

2) Envelope lines are drawn along the three sides of the triangular 
data space. The soil line is drawn as a straight line emanating from 
the origin. (For the assumption of no shadows and constant vegetation 
reflectivities, all sides of the triangle must be drawn straight.) 

3) Along the soil line, m is assumed equal to zero. Likewise, at 
the top of the triangle, m is assumed equal to one. 

4) The calculation of m for any data point within the triangle 
depends on the assumptions made. For constant vegetation reflectivity, 

m is assumed linearly related to the distance between the top and base of 
the triangle. For example, for a point exactly halfway between the top 
and the base, m is assumed equal to 50% cover. For the case of vegeta- 
tion reflectivity linearly related to percent cover, m is estimated on 
the basis of the shape of the envelope curves, which can be approximated 
as second order polynomials in m. 

The graphical results of the analyses applied to the Taos data are 
shown in Exhibits 20 through 23. Exhibits 20 and 21 include identical 
plots of all the data in the particular segment chosen (2250 pixels). It 
is noted that the data space takes on a classical triangular shape. 

Exhibit 20 also includes straight line envelope curves under the assump- 
tion of constant vegetation reflectivity. Exhibit 21 includes curved 
envelope lines under the assumption that vegetation reflectivity is 
linearly related to m. Percent cover for each case was then estimated 
according to the procedure outlined above. 

Exhibits 22 and 23 include only the data points with ground truth. 

In order to test the validity of the method, it was necessary to locate / 
the data points with ground truth within the triangle and compare the 
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actual percent cover, m, to that estimated graphically, m . 

& 

values resulting from the regression of m with m are shown 

S 

Table 5 

Results of Two-Band Models 

Two-Band Model Version R 2 

1) Constant Vegetation Reflectivity 0.34 

2) Variable r 

v 


The R 2 
in Table 5. 


0.30 
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Conclusions 

There are several important conclusions which can be drawn from this 
past year's work. In general, the results have established the feasibil- 
ity of estimating percent canopy cover using either physically-based or 
empirical models. The most successful empirical model for the Taos Study 
Area was the normalized vegetation index, which explained about 60 
percent of its variability in terms of m. About 50 percent of the 
variability of the direct beam equation for the visible band could be 
explained in terms of m using constant canopy and soil reflectivities. 

For the near infrared band, only 10 percent of the variability could be 
explained for the same conditions. The two-band model was able to 
explain about 35 percent of its variability in terms of m. 

While the use of NVI produced the best correlation, the results have 
depended on the availability of ground truth. Its application to other 
natural regions is therefore unknown. That is also true for the 
physically-based models which use only one band, even with the 
incorporation of slope and shadows. 

It has also been shown that the characteristic triangular shape of 
the two-band data space (MSS bands 2 and 4) can be explained in terms of 
percent canopy cover, although other vegetation parameters, such as leaf 
area and canopy reflectance, can also contribute to such shapes. A 
simple version assuming constant canopy reflectance provided encouraging 
results when applied to the Taos Study Area. 

The major advantage of the two-band model is in its ability to 
estimate canopy cover without ground truth. It uses the entire data set 
and provides the possibility of estimating other vegetation parameters by 
careful incorporation of the limits of the data space into the analysis. 



The simple version presented can be extended relatively easily to include 
topography, shadows, leaf area, and variability (randomness) in canopy 
reflectance (see Future Work), 

There is always some uncertainty regarding the error introduced 
through the data themselves. That includes inaccuracies in the estima- 
tion of percent cover from the aerial photographs, errors in registration 
of both the Landsat scene as well as the aerial photographs, anomalies in 
the Landsat data, errors introduced through data filtering, and finally, 
in the conversion of DN values to actual radiometric quantities. There 
appears to be no easy way to estimate the magnitude of each of these 
effects. Several of the more successful procedures should be applied to 
an entirely new data set outside the Taos Study Area. 
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Future Work 

There are at' least two research areas which have not yet been 
exploited in the development of a working canopy cover algorithm using 
Landsat data. The first involves an investigation of the random 
distribution of the radiances when plotted in two-, three-, or four-band 
space. The second involves the incorporation of additional surface 
parameters into the canopy cover algorithm. Future work will try to 
combine both these areas in order to extend the two-band procedure 
introduced this past year. Further details are described below together 
with an outline of the proposed mathematical formulation. 

While the two-band method incorporates the limits of the data space 
in determining m, it does not take full advantage of the distribution of 
the radiances within the data space. A hypothetical example is shown in 
Exhibit 24. By treating the observed radiances of the visible and near 
infrared bands (R2 ,R4) as random fields, a joint probability distribu- 
tion can be obtained. Most of the information of the segment is located 
within the interior of the data space, and not along the limits which 
envelope the data. Up to now, the simple two-band model has incorporated 
only information about the limits in a deterministic sense. By incorp- 
orating randomness into the formulation, information contained in the 
entire data set can be used. The above approach can be easily extended 
to include the randomness of three- and four-band data spaces. 

The second area which needs further research is the incorporation of 
additional surface features into the canopy cover algorithm. Much of the 
Taos Study Area is located in rolling foothills. It is obvious from the 
aerial photographs that vegetation type and percent cover is a function 
of the elevation, slope and azimuth of the land. The combined effect of 
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all those features is qualitatively seen in the Landsat images of 
individual bands. It has been shown earlier in this report that such 
parameters as shadows, ground slope and azimuth, leaf area index, and 
species can be incorporated into the physical model. Several of those 
features which were examined individually now need to be incorporated 
into one model. A specific problem which needs to be worked out is the 
effect of shadows for sloping surfaces on observed reflectivity. 

The following mathematical formulation is intended to incorporate 
the effects of both randomness and additional surface parameters into a 
physically based canopy cover algorithm. The formulation is presented 
for the case of two bands. The starting point is the direct beam 
equation which can be written very generally in terms of reflectances for 
the visible (MSS band 2) and near infrared (MSS band 4) wavelengths, 

r 2 (m,p,L 2 ,g,e,s,u,w) = m r 2v (p,L 2 , g ,0 ,s ,u) + (l-m)r 2s (g ,0 ,n,w) (23) 

r 4 (m,p,L 4> g,e,s,u,w) - m r 4v (p>L 4 ,g,0 ,s,u) + (l-m)r 4s (g ,0 ,n,w) (24) 

where 

m = percent vegetation cover 
p = percent shadow parameter 

= leaf area indices for bands 2 and 4, respectively 
g * ground slope 

0 = zenith angle, which includes effects of ground slope, ground 
azimuth, (altitude, hour angle, etc. ) 
s = vegetation species 
n = soil moisture 

u * a general parameter including factors other than those mentioned 
above causing variability in canopy reflectance, for example. 


stress, disease, etc 


30 


w = a general parameter including factors other than those listed 
above causing soil reflectance variability, including organic 
content, grain size and distribution, mineral content, ... etc. 

Written in the above form, all the variability in r^ and r^ 
is attributed to the variability in r^ v and For example, 

since it has been demonstrated earlier that up to 60 percent of the 
variability in r 2 can be explained by m, it follows that the balance 
or 40 percent is due to r 2v (p, 0, 0, s, u) and r 2g (0, 0, r, w) 

variability. 

The next step is to try to remove some of that variability in r 2v 
and r 2g by incorporating the effects of additional parameters, for 
instance, slope and zenith angle. Introducing a new function, g(g,e) to 
account for those effects, it can be shown that for flat surfaces (not 
necessarily horizontal), 

g(0>6) = cose/cosg 

where cos0 is given by equation (22). Therefore, 

r 2 = g(0,0)[m r' v (p,L 2 ,s,u) + (l-m)r£ g (n,w)] (26) 



■ g(6,0)[m r^ v (p,L 4 ,s,u) + (l-m)r^ g (r,w) 

(27) 

where 

r 2v " r 2 v /g(e, e) 

(28) 


r 2s = r 2s /g(fJ ’ 0) 

(29) 


r 4v = r 4v /g(3 ’ 0) 

(30) 


r 4s = r rs /g(6 » 0) 

(31) 
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Similarly, variability due to shadows can be removed from r^ v and 
r^ v using various geometric models indicated earlier. Let h(p ) represent 
the effect of shadows manifested by a geometric model on a sloping bed. 
The variable p is a function of geometric shape (i.e., cone, cylinder) as 
well as slope, azimuth and zenith angle. The reflectances can be 
rewritten 


r 2 “ g(8 >9) [m h(p) r" v (L 2 s,u) + (1-m) r£ g (n,w)] 
“ g(B,9)[m h( p ) r" v (L 4 s,u) + (1-m) r^ g (n,w)] 

where 

r 2v = r 2v /[g(6,0) 

r 4v “ r 4 v /tg(B, 0 ) h( p )] 


(32) 

(33) 


(34) 


In theory, one can continue developing functions f(L), j(s) and k(n) 
to remove more and more of the variability from ^ , r 4v* r 2s * anc * r 4s # 

At each stage, however, additional unknowns are introduced into the 
formulation, each associated with more ground truth. Since the intention 
is to develop a model without the need for ground truth, it is reasonable 
to stop at this point where no field data are required and use equations 
(32) and (33) to estimate m in conjunction with the available informa- 
tion. Let us first summarize the data and the major assumptions. 

The data include r^ and r^ obtained from the Landsat digital 
counts. It is assumed that the atmosphere is optically thin, thereby 
neglecting the scattering terms in the atmospheric radiative transfer 
equation. The parameters g and 0 and hence g(g,0) can be estimated from 
knowing time of Landsat overpass and from examination of USGS topographic 


maps 



3 


n 

l. 


The function h(p) can be obtained analytically and requires 
additional reasonable assumptions. Those include the shape of the 
canopy (i.e., cones) and the spatial distribution of trees. For natural 
conifer forests, a poisson distribution has often been assumed (Matern, 
1960, Strahler et al., 1981, Li et al., 1985). 

The unknowns in equations (32) and (33) are thus m, r”^, r^, r 4v > 
and r" . Making further the assumption 


r" v = f(x 2 ) A(s,u) (35) 

r" v = f(x 4 ) A(s,u) (36) 

r 2s “ g(X 2 ) B(n>w) (37) 

r As = g(X 4 ) B(n » w) (38) 


where f(X) and g(X) are normalized reflectivity functions ranging from 
zero to one, equations (32) and (33) become 

r 2 = g(8>0)[m h(p) f(X 2 ) A(s,u) + (1-m) g(X 2 ) B(n,w)] (39) 

r 4 = g(8,6)[m h( p ) f(X 4 ) A(s,u) + (1-m) g(X 4 ) B(n,w)] (40) 

There are now three unknowns, m, A(s,u) and B(n,w), associated with 
equations (39) and (40). While the problem still requires one more equa- 
tion, it has been greatly simplified by making only a few reasonable 
assumptions and without ground truth, except for topographic data. 

There are several approaches one can take to solve equations (39) 
and (40). The first is to obtain a derived probability distribution of m 
of the form 

m = mtr 2 , r 4 , A(s,w), B(n,w)] 


(41) 
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This can be done numerically. The random functions A(s,w) and B(n,w) 
might be obtained by conditionalizing on r^ and r^. 

A second approach is to introduce a third equation r^ or r^ using 
either of the two remaining Landsat MSS bands, although problems of 
independence between bands 1 and 2 or 3 and 4 might arise. 

A third approach is to simply assume a priori that the data falling 
on the lower base of the data space triangle reflect background. This 
leads to a direct calculation of the probability density function for 
B(n,w), by conditionalizing along that background line. Once determined, 
one can compute by means of equations (39) and (40), the two remaining 
unknowns, A(s,u) and m. 

It is not known which of the above three approaches work best for 
obtaining the desired third equation. The third approach, assuming a 
soil line, appears the most straightforward. It is also consistent with 
other research which attempts to use the soil line as an important 
parameter in vegetation modeling. Future work will examine and evaluate 
each of those alternatives. 
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